14. 重回帰分析
https://gyazo.com/4f7dfd91756d579cff6ea820ee6ff329
14.1. 多変量データ
世界47カ国の1961年から1965年における食物供給量と1974年の大腸がんの訂正死亡率のデータ
https://gyazo.com/dca75b677b2d539f3e6a5f2e159d2c6d
「総熱量」「肉類」「乳製品」「酒類」は1人1日あたりの熱量(kcal)
「大腸がん」は年齢構成によって訂正した人口10万人あたりの大腸がんによる死亡者数
以後「大腸がん」データと呼ぶ
ここでは1つの国に対して5種類の測定をしている
1つの観測対象から3回以上の測定を行ったデータ
14.1.1. 多変量散布図
対応あるデータの場合は散布図を描くことによってデータの様子を視覚的に調べることができた
多変量データの場合は多変量散布図を描くことによって、データの様子を視覚的に調べることができる https://gyazo.com/7bdcd2c9ddff801de0eb06898855fd9d
一般的に何らかの要素を四角形に並べたもの
左下から右下に向かっての斜めの要素は行番号と列番号が同じ
それ以外の行番号と列番号が異なる要素
$ i行$ j列の要素と$ j行$ i列の要素の内容が同じである行列
対称行列の形式で非対角要素に散布図を配したグラフ
対角要素には変数名(場合によってはヒストグラム)を置く
散布図を観察
「大腸がん」ともっとも正の相関関係が強いのは「肉類」であることがわかる
食品間の相関では「総熱量」と「肉類」の正の相関関係が強い
これはそもそも肉類の熱量が大きいことが原因と言えよう
14.1.2. 共分散行列・相関行列
多変量データに含まれる変数$ iと変数$ jの共分散を非対角要素の$ i行$ j列に並べた対角行列
共分散行列の対角要素には分散を配する
table: 表14-2 「大腸がんデータ」の共分散行列
総熱量 肉類 乳製品 酒類 大腸がん
総熱量 250240 74906 57018 23564 1672
肉類 74906 37514 22301 6497 746
乳製品 57018 22301 25376 2747 417
酒類 23564 6497 2747 5518 197
大腸がん 1672 746 417 197 20.43
多変量データに含まれる変数$ iと変数$ jの相関係数を非対角要素に並べた対称行列
相関行列の対角要素は、変数自身の相関だから$ 1.00
table: 表14-3 「大腸がんデータ」の平均・標準偏差・相関行列
総熱量 肉類 乳製品 酒類 大腸がん
平均 2871 307 243 109 8.54
標準偏差 500 194 159 74.3 4.52
総熱量 1.000 0.773 0.716 0.634 0.739
肉類 0.773 1.000 0.723 0.452 0.853
乳製品 0.716 0.723 1.000 0.232 0.579
酒類 0.634 0.452 0.232 1.000 0.588
大腸がん 0.739 0.853 0.579 0.588 1.000
「大腸がん」と最も正の相関関係が強いのは「肉類」($ 0.853)であることはすでに散布図から観察されていた
2番目に「大腸がん」と正の相関関係が強いのは「総熱量」($ 0.739)であることがわかる
14.2. 重回帰モデル
「大腸がん」を基準変数$ yとし、その死亡者の高低の予測を考える
「総熱量」「肉類」「乳製品」「酒類」のすべてを予測変数とする($ x_1から$ x_4)
(13.1)式に相当する関数としては、予測変数の重み付き和を選び、$ i番目の国の大腸がんを次のように予測する
$ \hat y_i = a + b_1x_{i1} + b_2x_{i2} + b_3x_{i3} + b_4x_{i4} \qquad (14.1)
複数の予測変数のによる基準変数の予測式
$ \hat y_i = a + b_1x_{i1} + \cdots + b_jx_{ij} + \cdots + b_px_{ip} \qquad (14.2)
添え字$ jは予測変数を示す
予測変数の数は$ p個
(14.2)式を利用して、$ y_iと$ x_{ij}の関係を分析する方法
重回帰式には基準変数$ y_iが登場しない
予測値が測定値にピタリと一致することは期待できないから
そこで(13.3)式に相当する
$ y_i = a + b_1x_{i1} + \cdots b_jx_{ij} + \cdots + b_px_{ip} + e_i \qquad (14.3)
を利用する
誤差変数$ e_iが、平均$ 0、標準偏差$ \sigma_eの正規分布に従い、$ e_iと$ x_{ij}が独立であるとすると、$ y_iの分布は正規分布の密度関数を利用して次のように表現される
$ f(y_i|\bm\theta) = f(y_i|a + b_1x_{i1} + \cdots + b_jx_{ij} + \cdots + b_px_{ip}, \sigma_e), \qquad (14.4)
$ ただし \bm\theta = (a, b_1, \cdots, b_j, \cdots, b_p, \sigma_e)
尤度は(13.8)式、事後分布は(13.11)式となり、MCMCを利用することにより、母数の爾後分布・生成量の事後分布・予測分布に従う乱数を生成することが可能になる
予測値の事後分布・重回帰式の事後分布・事後予測分布は、それぞれ(13.12)式・(13.13)式・(13.17)式のなかの$ a^{(t)} + b^{(t)}x_iを
$ a^{(t)} + b_1^{(t)}x_{i1} + \cdots + b_j^{(t)}x_{ij} + \cdots + b_p^{(t)}x_{ip} \qquad (14.5)
に置き換えればよい
決定係数の事後分布は(13.16)式のまま
14.2.1. 標準偏回帰係数
単回帰分析とは異なり、予測変数が複数ある重回帰分析では重みの大きさを互いに比較する必要が生じる
しかし偏回帰係数の大きさは、単純に比較することはできない
たとえば$ x_1「総熱量」の標準偏差は$ 500であり、$ x_4「酒類」の標準偏差は$ 74.3
仮に2つの変数の偏回帰係数の値が同じであったとしても、変数の散らばりの大きさが違うから、基準変数に対する影響力が異なり、比較は困難
基準変数と予測変数を、平均$ 0、標準偏差$ 1の標準得点に変換した時の偏回帰係数
標準得点に変換して標準偏回帰係数を求める
$ y_iあるいは$ x_{ij}を含まない項を$ \cdotsで表現すると(14.3)式は以下のように表現できる
$ \begin{aligned} y_i & = b_jx_{ij} + \cdots \\ & \small{[両辺に母平均などで同じ項を足して引き]} \\ y_i - \mu_y + \mu_y & = b_jx_{ij} - b_j\mu_{x_j} + b_j\mu_{x_j} + \cdots \\ & \small{[左辺 +\mu_y を移項し、右辺 +b_j\mu_{x_j}とともに\cdotsに含め]} \\ y_i - \mu_y & = b_jx_{ij} - b_j\mu_{x_j} + \cdots \\ & \small{[両辺をそれぞれの標準偏差で割り、そして掛け]} \\ \frac{y_i - \mu_y}{s_y}s_y & = b_js_{x_j}\frac{x_{ij}-\mu_{x_j}}{s_{x_j}} + \cdots \\ & \small{[平均を引き標準偏差で割った部分を標準得点で置き換え]} \\ z_{yi}s_y & = b_js_{x_j}z_{xij} + \cdots \\ & \small{[両辺をs_yで割ると]} \\ z_{yi} & = \frac{b_js_{x_j}}{s_y}z_{xij} + \cdots \end{aligned}
以上のことから標準偏回帰係数$ b_j^*は以下のように導かれる
$ b_j^* = \frac{b_js_{x_j}}{s_y}, \quad (j=1, \cdots, p) \qquad (14.6)
切片$ aは、値が$ 1だけの予測変数にかかる重みと解釈することができる
この場合標準偏差は$ 0になるから, (14.6)式の分子が$ 0となる
したがって標準得点の重回帰式では切片は常に$ 0となる
標準回帰係数の事後分布のための生成量は(13.14)式を考慮し、
$ b_j^{*(t)} = \frac{b_j^{(t)}s_{x_j}}{s_y} \qquad (14.7)
によって近似する
14.2.2. 重相関係数
予測変数による基準変数の予測の精度としては、すでに決定係数を学んだ
基準変数の測定値と予測値の相関係数
$ y = \hat y + eであることや$ \hat yと$ eが互いに独立であることや、(7.17)式の関係等を利用すると、重相関係数$ p_{y\hat y}は
$ \rho_{y \hat y} = \frac{\sigma_{y \hat y}}{\sigma_y \sigma_{\hat y}} = \frac{cov(\hat y + e, \hat y)}{\sigma_y\sigma_{\hat y}} = \frac{\sigma_{\hat y}^2}{\sigma_y\sigma_{\hat y}} = \frac{\sigma_{\hat y}}{\sigma_y} = \eta \qquad (14.8)
のように決定計数(13.15)式の平方根に一致する
ここで関数$ cov(\quad)は共分散を計算する関数
したがって、重相関係数の事後分布は、生成量
$ \rho_{y \hat y}^{(t)} = \eta^{(t)} \qquad (14.9)
によって近似できる
14.3. 分析結果
「大腸がん」を基準変数とし、「総熱量」「肉類」「乳製品」「酒類」を予測変数(それぞれ$ x_1から$ x_4)とし、(14.1)式で重回帰分析を行う
$ 21000個の乱数を5本発生させ、バーンイン期間を$ 1000とし、$ T=100000の乱数によって母数の事後分布を近似した
14.3.1. 変数選択
母数の推定結果
table: 表14-4 重回帰モデルの母数と決定係数と重相関係数の事後分布
EAP post.sd 2.5% 5% 50% 95% 97.5%
a 0.58332 3.00036 -5.32721 -4.34798 0.58711 5.50150 6.48870
b_1 0.00052 0.00143 -0.00230 -0.00182 0.00051 0.00287 0.00334
b_2 0.01732 0.00299 0.01146 0.01243 0.01731 0.02227 0.02324
b_3 -0.00148 0.00361 -0.00864 -0.00744 -0.00146 0.00445 0.00563
b_4 0.01394 0.00648 0.00120 0.00327 0.01397 0.02444 0.02659
σ_e 2.31176 0.26156 1.86523 1.92531 2.28817 2.77714 2.88811
4つの偏回帰係数のなかで最初に注意すべきは「『乳製品』の偏回帰係数$ b_3のEAPが$ -0.00148であり、負である」こと
ここだけを見て「『乳製品』の増加は『大腸がん』の減少に寄与する」と解釈してはいけない
図14-1や表14-3を観察すれば明らかなように、「乳製品」と「大腸がん」には正の相関関係がある
偏回帰係数は「自身以外の予測変数が一定であるときに、当該予測変数が1単位動いたときの、基準変数の平均的変化量」と解釈する
「乳製品」の$ b_3のEAPが$ -0.00148であるということは「『総熱量』『肉類』『酒類』が一定であるときに、『乳製品』を$ 100kcal増加させると『大腸がん』は$ 0.148(= 0.00148 \times 100)人減少する」と解釈するのが正しい
しかし、平均が$ 243kcalで標準偏差が$ 159kcalの「乳製品」をさらに$ 100kcal増加させることは大変なこと
対して減るのはたった$ 0.15人
乳製品の偏回帰係数$ b_3には実質的な効果を期待できない程度といえる
偏回帰係数の正しい解釈は難しい
例えば「総熱量」を一定にして「乳製品」を$ 100kcal増加させることは、想像しにくい
予測変数の数$ pが大きくなると「自身以外の予測変数が一定であるときに」という前提が実質科学的に意味をもたなくなるケースが生じやすくなる
そうなったら偏回帰係数の解釈はあきらめ、分析の主目的を基準変数の予測とする
$ pが大きい場合は、小さい場合に比べて決定係数や重相関係数の点推定値が大きくなり、予測精度が向上しているように見えることがある
簡便なためにしばしば利用される確信区間を利用した方法
$ b_2は95%の確信で「$ 0.01243以上である」といえる
$ b_4は95%の確信で「$ 0.00327以上である」といえる
この2つの母数は、正の値である確率が高い
そこで「総熱量」「乳製品」を削除し、「肉類」「酒類」を残すという変数選択を行う
予測変数を「肉類」「酒類」とし、$ p=2で重回帰分析を再び行い、母数と決定係数と重相関係数の推定をだす
table: 表14-5 変数選択後の母数と決定係数と重相関係数の事後分布
EAP post.sd 2.5% 5% 50% 95% 97.5%
a 1.571 0.681 0.240 0.459 1.569 2.695 2.922
b_2 0.017 0.002 0.013 0.014 0.017 0.020 0.021
b_4 0.016 0.005 0.006 0.007 0.016 0.024 0.026
σ_e 2.261 0.250 1.837 1.894 2.239 2.708 2.817
決定係数 0.756 0.051 0.637 0.662 0.762 0.827 0.837
重相関係数 0.869 0.030 0.798 0.813 0.873 0.910 0.915
b_2^* 0.737 0.082 0.574 0.602 0.737 0.872 0.899
b_4^* 0.256 0.083 0.093 0.119 0.256 0.392 0.419
表14-4のWAICは$ 215.6であり、表14-5のWAICは$ 211.8で、WAICの観点からも変数選択は支持されている
14.3.2. 偏回帰係数・決定係数・重相関係数
切片の推定値は$ 1.571(0.681)[0.240, 2.922] であり、「肉類」と「酒類」が$ 0であるときの「大腸がん」の予測値
菜食主義でアルコールを嗜まない場合の死亡者数は約1.6人と予測される
$ b_2 の推定値は$ 0.017(0.002)[0.013, 0.021] であり、「酒類」が一定であるときに、「肉類」が$ 100kcal増加すると死亡者が平均的に約$ 1.7人増加する
$ b_4 の推定値は$ 0.016(0.005)[0.006, 0.026] であり、「肉類」が一定であるときに「酒類」が$ 30kcal増加すると、死亡者が平均的に約$ 0.48人増加する
決定係数の推定値は$ 0.756(0.051)[0.637, 0.837] であり、基準変数の変動は予測変数の組によって75.6%説明された
重相関係数の推定値は$ 0.869(0.030)[0.798, 0.915] であり、基準変数の測定値と予測値の相関は約$ 0.87
$ b_2^* の推定値は$ 0.737(0.082)[0.574, 0.899] であり、$ b_4^* の推定値は$ 0.256(0.083)[0.093, 0.419]
このことから「酒類」より「肉類」の方が標準偏回帰係数が大きいと解釈する
$ b_2^*と$ b_4^*の差の生成量を構成し、その事後分布から差に関する推測を行うことも可能
14.3.3. 回帰平面
3次元空間に3つの変数の測定点をプロットしたグラフ
https://gyazo.com/d32e55fbbe7f89e96203370c310db00b
立体的なイメージを促すために測定点から$ xy平面に垂線を下している
$ p=2の重回帰式は3次元空間内に1枚の平面を張る
図14-2には回帰平面
$ 「大腸がん」= 1.571 + 0.017「肉類」 + 0.016「酒類」 \qquad (14.10)
も描いている。
回帰平面の近くの上下に測定点がまとわりついている
14.3.4. 残差プロット
単回帰分析では予測変数×残差の残差プロットを描くことが有効であった
重回帰分析では図14-3のように、基準変数の予測値×残差の残差プロットを描くことが有効
https://gyazo.com/1a03fc316966f8d70544375c923a3a6f
たとえば番号5の国が予測値より死亡者が$ 6人多いこと、番号43の国が予測値より死亡者数が約$ 4人少ないことが観察される
14.3.5. 予測分布
単回帰分析では、広い範囲の間隔の短い等差数列を$ x^*に与えることにより事後予測分布を近似した
しかし、重回帰分析では予測変数の次元数が高いので網羅的に予測変数の状態を準備するのは手間がかかる
そこで予測変数のいくつかの測定点に対する基準変数の事後予測分布を考察することが多い
table: 表14-6 予測値・事後標準偏差・事後分布95%点・事後予測分布95%点
肉類 酒類 予測値 post.sd 事後分布95%点 sd 事後予測分布95%点
0.0 0.0 1.6 0.7 2.7 2.4 5.4
300 0.0 6.7 0.6 7.8 2.4 10.6
0.0 100 3.1 0.7 4.2 2.4 7.0
150 50 4.9 0.5 5.7 2.3 8.7
たとえば$ x_2 = 150, x_4 = 50における$ y^*の死亡者の予測値は$ 4.9人であり、上側点は$ 8.7人である